/*
Adapted from AAPlus library by PJ Naughter
Permission was granted to convert the library to C#

This source is modified from PJ Naugters AA+ lib
(http://www.naughter.com/aa.html)
and therefore is covered under his copyright license, which
is as follows:

Purpose: Implementation for the algorithms for an elliptical orbit
Created: PJN / 29-12-2003
History: PJN / 24-05-2004 1. Fixed a missing break statement in Calculate. Thanks to
                          Carsten A. Arnholm for reporting this bug. 
                          2. Also fixed an issue with the calculation of the apparent distance to 
                          the Sun.
         PJN / 31-12-2004 1. Fix for MinorPlanetMagnitude where the phase angle was
                          being incorrectly converted from Radians to Degress when it was already
                          in degrees. Thanks to Martin Burri for reporting this problem.
         PJN / 05-06-2006 1. Fixed a bug in Calculate(double JD, EllipticalObject object)
                          where the correction for nutation was incorrectly using the Mean obliquity of
                          the ecliptic instead of the true value. The results from the test program now 
                          agree much more closely with the example Meeus provides which is the position 
                          of Venus on 1992 Dec. 20 at 0h Dynamical Time. I've also checked the positions
                          against the JPL Horizons web site and the agreement is much better. Because the
                          True obliquity of the Ecliptic is defined as the mean obliquity of the ecliptic
                          plus the nutation in obliquity, it is relatively easy to determine the magnitude
                          of error this was causing. From the chapter on Nutation in the book, and 
                          specifically the table which gives the cosine coefficients for nutation in 
                          obliquity you can see that the absolute worst case error would be the sum of the 
                          absolute values of all of the coefficients and would have been c. 10 arc seconds 
                          of degree, which is not a small amount!. This value would be an absolute worst 
                          case and I would expect the average error value to be much much smaller 
                          (probably much less than an arc second). Anyway the bug has now been fixed. 
                          Thanks to Patrick Wong for pointing out this rather significant bug. 
         PJN / 10-11-2008 1. Fixed a bug in Calculate(double JD, 
                          const CAAEllipticalObjectElements& elements) in the calculation of the 
                          heliocentric rectangular ecliptical, the heliocentric ecliptical latitude and 
                          the heliocentric ecliptical longitude coordinates. The code incorrectly used the 
                          value "omega" instead of "w" in its calculation of the value "u". Unfortunately 
                          there is no worked examples in Jean Meeus's book for these particular values, 
                          hence resulting in my coding errors. Thanks to Carsten A. Arnholm for reporting 
                          this bug. 

Copyright (c) 2003 - 2009 by PJ Naughter (Web: www.naughter.com, Email: pjna@naughter.com)

All rights reserved.

Copyright / Usage Details:

You are allowed to include the source code in any product (commercial, shareware, freeware or otherwise) 
when your product is released in binary form. You are allowed to modify the source code in any way you want 
except you cannot modify the copyright details at the top of each module. If you want to distribute source 
code with your application, then you are only allowed to distribute versions released by the author. This is 
to maintain a single distribution point for the source code. 

*/

using System;
using System.Diagnostics;

namespace AAPlus
{
    public class EllipticalPlanetaryDetails
    {
        public double ApparentGeocentricLongitude = 0.0;
        public double ApparentGeocentricLatitude = 0.0;
        public double ApparentGeocentricDistance = 0.0;
        public double ApparentLightTime = 0.0;
        public double ApparentGeocentricRA = 0.0;
        public double ApparentGeocentricDeclination = 0.0;
    };

    public class EllipticalObjectElements
    {
        public double a = 0.0;
        public double e = 0.0;
        public double i = 0.0;
        public double w = 0.0;
        public double omega = 0.0;
        public double JDEquinox = 0.0;
        public double T = 0.0;
    };

    public class EllipticalObjectDetails
    {
        public Coordinate3D HeliocentricRectangularEquatorial;
        public Coordinate3D HeliocentricRectangularEcliptical;
        public double HeliocentricEclipticLongitude = 0.0;
        public double HeliocentricEclipticLatitude = 0.0;
        public double TrueGeocentricRA = 0.0;
        public double TrueGeocentricDeclination = 0.0;
        public double TrueGeocentricDistance = 0.0;
        public double TrueGeocentricLightTime = 0.0;
        public double AstrometricGeocenticRA = 0.0;
        public double AstrometricGeocentricDeclination = 0.0;
        public double AstrometricGeocentricDistance = 0.0;
        public double AstrometricGeocentricLightTime = 0.0;
        public double Elongation = 0.0;
        public double PhaseAngle = 0.0;

        public EllipticalObjectDetails()
        {
            HeliocentricRectangularEquatorial = new Coordinate3D();
            HeliocentricRectangularEcliptical = new Coordinate3D();
        }
    };

	public static class Elliptical
	{
		public enum EllipticalObject
		{
			Sun,
			Mercury,
			Venus,
			Mars,
			Jupiter,
			Saturn,
			Uranus,
			Neptune,
			Pluto
		};
		
public static double DistanceToLightTime(double Distance)
{
  return Distance * 0.0057755183;
}

public static EllipticalPlanetaryDetails Calculate(double JD, EllipticalObject obj)
{
  //What will the the return value
  EllipticalPlanetaryDetails details = new EllipticalPlanetaryDetails();

  double JD0 = JD;
  double L0 = 0;
  double B0 = 0;
  double R0 = 0;
  double cosB0 = 0;
  if (obj != EllipticalObject.Sun)
  {
    L0 = Earth.EclipticLongitude(JD0);
    B0 = Earth.EclipticLatitude(JD0);
    R0 = Earth.RadiusVector(JD0);
    L0 = CoordinateTransformation.DegreesToRadians(L0);
    B0 = CoordinateTransformation.DegreesToRadians(B0);
    cosB0 = Math.Cos(B0);
  }


  //Calculate the initial values
  double L = 0;
  double B = 0;
  double R = 0;
  switch (obj)
  {
    case EllipticalObject.Sun:
    {
      L = Sun.GeometricEclipticLongitude(JD0);
      B = Sun.GeometricEclipticLatitude(JD0);
      R = Earth.RadiusVector(JD0);
      break;
    }
    case EllipticalObject.Mercury:
    {
      L = Mercury.EclipticLongitude(JD0);
      B = Mercury.EclipticLatitude(JD0);
      R = Mercury.RadiusVector(JD0);
      break;
    }
    case EllipticalObject.Venus:
    {
      L = Venus.EclipticLongitude(JD0);
      B = Venus.EclipticLatitude(JD0);
      R = Venus.RadiusVector(JD0);
      break;
    }
    case EllipticalObject.Mars:
    {
      L = Mars.EclipticLongitude(JD0);
      B = Mars.EclipticLatitude(JD0);
      R = Mars.RadiusVector(JD0);
      break;
    }
    case EllipticalObject.Jupiter:
    {
      L = Jupiter.EclipticLongitude(JD0);
      B = Jupiter.EclipticLatitude(JD0);
      R = Jupiter.RadiusVector(JD0);
      break;
    }
    case EllipticalObject.Saturn:
    {
      L = Saturn.EclipticLongitude(JD0);
      B = Saturn.EclipticLatitude(JD0);
      R = Saturn.RadiusVector(JD0);
      break;
    }
    case EllipticalObject.Uranus:
    {
      L = Uranus.EclipticLongitude(JD0);
      B = Uranus.EclipticLatitude(JD0);
      R = Uranus.RadiusVector(JD0);
      break;
    }
    case EllipticalObject.Neptune:
    {
      L = Neptune.EclipticLongitude(JD0);
      B = Neptune.EclipticLatitude(JD0);
      R = Neptune.RadiusVector(JD0);
      break;
    }
    case EllipticalObject.Pluto:
    {
      L = Pluto.EclipticLongitude(JD0);
      B = Pluto.EclipticLatitude(JD0);
      R = Pluto.RadiusVector(JD0);
      break;
    }
    default:
    {
      Debug.Assert(false);
      break;
    }
  }

  bool bRecalc = true;
  bool bFirstRecalc = true;
  double LPrevious = 0;
  double BPrevious = 0;
  double RPrevious = 0;
  while (bRecalc)
  {
    switch (obj)
    {
      case EllipticalObject.Sun:
      {
        L = Sun.GeometricEclipticLongitude(JD0);
        B = Sun.GeometricEclipticLatitude(JD0);
        R = Earth.RadiusVector(JD0);
        break;
      }
      case EllipticalObject.Mercury:
      {
        L = Mercury.EclipticLongitude(JD0);
        B = Mercury.EclipticLatitude(JD0);
        R = Mercury.RadiusVector(JD0);
        break;
      }
      case EllipticalObject.Venus:
      {
        L = Venus.EclipticLongitude(JD0);
        B = Venus.EclipticLatitude(JD0);
        R = Venus.RadiusVector(JD0);
        break;
      }
      case EllipticalObject.Mars:
      {
        L = Mars.EclipticLongitude(JD0);
        B = Mars.EclipticLatitude(JD0);
        R = Mars.RadiusVector(JD0);
        break;
      }
      case EllipticalObject.Jupiter:
      {
        L = Jupiter.EclipticLongitude(JD0);
        B = Jupiter.EclipticLatitude(JD0);
        R = Jupiter.RadiusVector(JD0);
        break;
      }
      case EllipticalObject.Saturn:
      {
        L = Saturn.EclipticLongitude(JD0);
        B = Saturn.EclipticLatitude(JD0);
        R = Saturn.RadiusVector(JD0);
        break;
      }
      case EllipticalObject.Uranus:
      {
        L = Uranus.EclipticLongitude(JD0);
        B = Uranus.EclipticLatitude(JD0);
        R = Uranus.RadiusVector(JD0);
        break;
      }
      case EllipticalObject.Neptune:
      {
        L = Neptune.EclipticLongitude(JD0);
        B = Neptune.EclipticLatitude(JD0);
        R = Neptune.RadiusVector(JD0);
        break;
      }
      case EllipticalObject.Pluto:
      {
        L = Pluto.EclipticLongitude(JD0);
        B = Pluto.EclipticLatitude(JD0);
        R = Pluto.RadiusVector(JD0);
        break;
      }
      default:
      {
        Debug.Assert(false);
        break;
      }
    }

    if (!bFirstRecalc)
    {
      bRecalc = ((Math.Abs(L - LPrevious) > 0.00001) || (Math.Abs(B - BPrevious) > 0.00001) || (Math.Abs(R - RPrevious) > 0.000001));
      LPrevious = L;
      BPrevious = B;
      RPrevious = R;
    }
    else
      bFirstRecalc = false;  

    //Calculate the new value
    if (bRecalc)
    {
      double distance = 0;
      if (obj != EllipticalObject.Sun)
      {
        double Lrad = CoordinateTransformation.DegreesToRadians(L);
        double Brad = CoordinateTransformation.DegreesToRadians(B);
        double cosB = Math.Cos(Brad);
        double cosL = Math.Cos(Lrad);
        double x = R * cosB * cosL - R0 * cosB0 * Math.Cos(L0);
        double y = R * cosB * Math.Sin(Lrad) - R0 * cosB0 * Math.Sin(L0);
        double z = R * Math.Sin(Brad) - R0 * Math.Sin(B0);
        distance = Math.Sqrt(x*x + y*y + z*z);
      }
      else
        distance = R; //Distance to the sun from the earth is in fact the radius vector

      //Prepare for the next loop around
      JD0 = JD - DistanceToLightTime(distance);
    }
  }

    {
        double Lrad = CoordinateTransformation.DegreesToRadians(L);
        double Brad = CoordinateTransformation.DegreesToRadians(B);
        double cosB = Math.Cos(Brad);
        double cosL = Math.Cos(Lrad);
        double x = R*cosB*cosL - R0*cosB0*Math.Cos(L0);
        double y = R*cosB*Math.Sin(Lrad) - R0*cosB0*Math.Sin(L0);
        double z = R*Math.Sin(Brad) - R0*Math.Sin(B0);
        double x2 = x*x;
        double y2 = y*y;

        details.ApparentGeocentricLatitude = CoordinateTransformation.RadiansToDegrees(Math.Atan2(z, Math.Sqrt(x2 + y2)));
        details.ApparentGeocentricDistance = Math.Sqrt(x2 + y2 + z*z);
        details.ApparentGeocentricLongitude =
            CoordinateTransformation.MapTo0To360Range(CoordinateTransformation.RadiansToDegrees(Math.Atan2(y, x)));
        details.ApparentLightTime = DistanceToLightTime(details.ApparentGeocentricDistance);

        //Adjust for Aberration
        Coordinate2D Aberration = AAPlus.Aberration.EclipticAberration(details.ApparentGeocentricLongitude,
                                                                       details.ApparentGeocentricLatitude, JD);
        details.ApparentGeocentricLongitude += Aberration.X;
        details.ApparentGeocentricLatitude += Aberration.Y;

        //convert to the FK5 system
        double DeltaLong = FK5.CorrectionInLongitude(details.ApparentGeocentricLongitude,
                                                     details.ApparentGeocentricLatitude, JD);
        details.ApparentGeocentricLatitude += FK5.CorrectionInLatitude(details.ApparentGeocentricLongitude, JD);
        details.ApparentGeocentricLongitude += DeltaLong;

        //Correct for nutation
        double NutationInLongitude = Nutation.NutationInLongitude(JD);
        double Epsilon = Nutation.TrueObliquityOfEcliptic(JD);
        details.ApparentGeocentricLongitude += CoordinateTransformation.DMSToDegrees(0, 0, NutationInLongitude);

        //Convert to RA and Dec
        Coordinate2D ApparentEqu = CoordinateTransformation.Ecliptic2Equatorial(details.ApparentGeocentricLongitude,
                                                                                   details.ApparentGeocentricLatitude,
                                                                                   Epsilon);
        details.ApparentGeocentricRA = ApparentEqu.X;
        details.ApparentGeocentricDeclination = ApparentEqu.Y;
    }
    return details;
}

public static double SemiMajorAxisFromPerihelionDistance(double q, double e)
{
  return q / (1 - e);
}

public static double MeanMotionFromSemiMajorAxis(double a)
{
  return 0.9856076686 / (a * Math.Sqrt(a));
}

public static EllipticalObjectDetails Calculate(double JD, out EllipticalObjectElements elements)
{
    elements = new EllipticalObjectElements();
    double Epsilon = Nutation.MeanObliquityOfEcliptic(elements.JDEquinox);

  double JD0 = JD;

  //What will be the return value
  var details = new EllipticalObjectDetails();

  Epsilon = CoordinateTransformation.DegreesToRadians(Epsilon);
  double omega = CoordinateTransformation.DegreesToRadians(elements.omega);
  double w = CoordinateTransformation.DegreesToRadians(elements.w);
  double i = CoordinateTransformation.DegreesToRadians(elements.i);

  double sinEpsilon = Math.Sin(Epsilon);
  double cosEpsilon = Math.Cos(Epsilon);
  double sinOmega = Math.Sin(omega);
  double cosOmega = Math.Cos(omega);
  double cosi = Math.Cos(i);
  double sini = Math.Sin(i);

  double F = cosOmega;
  double G = sinOmega * cosEpsilon;
  double H = sinOmega * sinEpsilon;
  double P = -sinOmega * cosi;
  double Q = cosOmega*cosi*cosEpsilon - sini*sinEpsilon;
  double R = cosOmega*cosi*sinEpsilon + sini*cosEpsilon;
  double a = Math.Sqrt(F*F + P*P);
  double b = Math.Sqrt(G*G + Q*Q);
  double c = Math.Sqrt(H*H + R*R);
  double A = Math.Atan2(F, P);
  double B = Math.Atan2(G, Q);
  double C = Math.Atan2(H, R);
  double n = MeanMotionFromSemiMajorAxis(elements.a);

  Coordinate3D SunCoord = Sun.EquatorialRectangularCoordinatesAnyEquinox(JD, elements.JDEquinox);

  for (int j=0; j<2; j++)
  {
    double M = n * (JD0 - elements.T);
    double E = Kepler.Calculate(M, elements.e);
    E = CoordinateTransformation.DegreesToRadians(E);
    double v = 2*Math.Atan(Math.Sqrt((1 + elements.e) / (1 - elements.e)) * Math.Tan(E/2));
    double r = elements.a * (1 - elements.e*Math.Cos(E));
    double x = r * a * Math.Sin(A + w + v);
    double y = r * b * Math.Sin(B + w + v);
    double z = r * c * Math.Sin(C + w + v);

    if (j == 0)
    {
      details.HeliocentricRectangularEquatorial.X = x;
      details.HeliocentricRectangularEquatorial.Y = y;
      details.HeliocentricRectangularEquatorial.Z = z;

      //Calculate the heliocentric ecliptic coordinates also
      double u = w + v;
      double cosu = Math.Cos(u);
      double sinu = Math.Sin(u);

      details.HeliocentricRectangularEcliptical.X = r * (cosOmega*cosu - sinOmega*sinu*cosi);
      details.HeliocentricRectangularEcliptical.Y = r * (sinOmega*cosu + cosOmega*sinu*cosi);
      details.HeliocentricRectangularEcliptical.Z = r*sini*sinu;

      details.HeliocentricEclipticLongitude = CoordinateTransformation.MapTo0To360Range(CoordinateTransformation.RadiansToDegrees(Math.Atan2(details.HeliocentricRectangularEcliptical.Y, details.HeliocentricRectangularEcliptical.X))); 
      details.HeliocentricEclipticLatitude = CoordinateTransformation.RadiansToDegrees(Math.Asin(details.HeliocentricRectangularEcliptical.Z / r));
    }

    double psi = SunCoord.X + x;
    double nu = SunCoord.Y + y;
    double sigma = SunCoord.Z + z;

    double Alpha = Math.Atan2(nu, psi);
    Alpha = CoordinateTransformation.RadiansToDegrees(Alpha);
    double Delta = Math.Atan2(sigma, Math.Sqrt(psi*psi + nu*nu));
    Delta = CoordinateTransformation.RadiansToDegrees(Delta);
    double Distance = Math.Sqrt(psi*psi + nu*nu + sigma*sigma);

    if (j == 0)
    {
      details.TrueGeocentricRA = CoordinateTransformation.MapTo0To24Range(Alpha / 15);
      details.TrueGeocentricDeclination = Delta;
      details.TrueGeocentricDistance = Distance;
      details.TrueGeocentricLightTime = DistanceToLightTime(Distance);
    }
    else
    {
      details.AstrometricGeocenticRA = CoordinateTransformation.MapTo0To24Range(Alpha / 15);
      details.AstrometricGeocentricDeclination = Delta;
      details.AstrometricGeocentricDistance = Distance;
      details.AstrometricGeocentricLightTime = DistanceToLightTime(Distance);

      double RES = Math.Sqrt(SunCoord.X*SunCoord.X + SunCoord.Y*SunCoord.Y + SunCoord.Z*SunCoord.Z);

      details.Elongation = Math.Acos((RES*RES + Distance*Distance - r*r) / (2 * RES * Distance));
      details.Elongation = CoordinateTransformation.RadiansToDegrees(details.Elongation);

      details.PhaseAngle = Math.Acos((r*r + Distance*Distance - RES*RES) / (2 * r * Distance));
      details.PhaseAngle = CoordinateTransformation.RadiansToDegrees(details.PhaseAngle);
    }

    if (j == 0) //Prepare for the next loop around
      JD0 = JD - details.TrueGeocentricLightTime;
  }
    
  return details;
}

public static double InstantaneousVelocity(double r, double a)
{
  return 42.1219 * Math.Sqrt((1/r) - (1/(2*a)));
}

public static double VelocityAtPerihelion(double e, double a)
{
  return 29.7847 / Math.Sqrt(a) * Math.Sqrt((1+e)/(1-e));
}

public static double VelocityAtAphelion(double e, double a)
{
  return 29.7847 / Math.Sqrt(a) * Math.Sqrt((1-e)/(1+e));
}

public static double LengthOfEllipse(double e, double a)
{
  double b = a * Math.Sqrt(1 - e*e);
  return CoordinateTransformation.PI() * (3 * (a+b) - Math.Sqrt((a+3*b)*(3*a + b)));
}

public static double CometMagnitude(double g, double delta, double k, double r)
{
  return g + 5*Math.Log10(delta) + k*Math.Log10(r);
}

public static double MinorPlanetMagnitude(double H, double delta, double G, double r, double PhaseAngle)
{
  //Convert from degrees to radians
  PhaseAngle = CoordinateTransformation.DegreesToRadians(PhaseAngle);

  double phi1 = Math.Exp(-3.33*Math.Pow(Math.Tan(PhaseAngle/2), 0.63));
  double phi2 = Math.Exp(-1.87*Math.Pow(Math.Tan(PhaseAngle/2), 1.22));

  return H + 5*Math.Log10(r*delta) - 2.5*Math.Log10((1 - G)*phi1 + G*phi2);
}

	}
}
